function sf = solve_forward(param,fp,fp_parfor,h)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Simulate the model;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

%Store model's parameters;
alpha = param(1:11);
TFP = param(12:13);   
gamma = param(17:21);

%Pre-define matrices to store simulated data;
sim_data = NaN( 3*fp_parfor.n_simulation , 13 ) ;
sim_data_network = NaN( 3*fp_parfor.n_simulation , 5 ) ;

grid =  cell(fp.periods);
grid_peers =  cell(fp.periods);
grid_inv0   = cell(fp.periods);
grid_inv1   = cell(fp.periods);
hc = cell((fp.periods-1),1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
% Draw Initial Conditions;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

init_draw_log_skills = fp_parfor.init_draws;           
hc{1}= exp(init_draw_log_skills);

Peers = NaN(size(hc{1},1),fp.periods);
distrib_skills   = NaN(size(hc{1},1),fp.periods);
distrib_PS = NaN(size(hc{1},1),fp.periods-1);
Skills = NaN(size(hc{1},1),fp.periods);
Peers_log = NaN(size(hc{1},1),fp.periods);

Skills(:,1) = hc{1};            
 
%Here loop over t for next periods ;
               
for t=1:(fp.periods-1)
                
        %Random graph generation;

        [H1, H2] = ndgrid(hc{t}, hc{t});          


        if t==1
        
        U_links1 = gamma(1).*ones(size(H1)) +  ...
                   gamma(2).*log(H1) + gamma(3).*log(H2) + gamma(4).*(log(H1) - log(H2)).^(2);                                                 
        U_links2 = gamma(1).*ones(size(H1)) +  ...
                   gamma(2).*log(H2) + gamma(3).*log(H1) + gamma(4).*(log(H1) - log(H2)).^(2) ;   
               
        else
            
        %t-1 PS stored in the loop;     
        [PS1, PS2] = ndgrid(PS, PS);              
            
        U_links1 = gamma(1).*ones(size(H1)) +  ...
                   gamma(2).*log(H1) + gamma(3).*log(H2) + gamma(4).*(log(H1) - log(H2)).^(2) + ...
                   gamma(5).*( (log(H1) - log(H2)).^(2) ).*(log(H1) - log(H2)>0).*PS1;                                                 
        U_links2 = gamma(1).*ones(size(H1)) +  ...
                   gamma(2).*log(H2) + gamma(3).*log(H1) + gamma(4).*(log(H1) - log(H2)).^(2) + ...
                   gamma(5).*( (log(H1) - log(H2)).^(2) ).*(log(H2) - log(H1)>=0).*PS2;               
            
        end
               
               
        Prob_Links = (1./(1+1./exp(U_links1))).*(1./(1+1./exp(U_links2))) ;   
                      
        ADJ = ( Prob_Links >= fp_parfor.peers_forward(:,:,1,t));        
        ADJ(logical(eye(size(ADJ)))) = 0;          

        peers_skills = hc{t} ;
         
        peers =  repmat(peers_skills',[size(Prob_Links,1) ,1]) ;      
        temp = ADJ.*peers; 
        H_peers =  nansum(temp,2)./sum(temp~=0 &  isnan(temp)==0,2) ;                                 
        Peers(:,t) = H_peers;                
        
        %Compute the mean-log peers for the indirect inference moments; 
        temp = ADJ.*log(peers); 
        H_peers_logs =  nansum(temp,2)./sum(temp~=0 &  isnan(temp)==0,2) ;                          
        Peers_log(:,t) = H_peers_logs;

        %Compute Statistics about the networks;

        %Number of Friends
        n_friends = sum( ADJ , 2 );

        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
        % Simulate Parental Optimal Behavior;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
                                               
        State_vars = [hc{t} ,  H_peers ];        
        value0 = Interp_function_simulation(fp_parfor.estim_value0{t},State_vars);
        value1 = Interp_function_simulation(fp_parfor.estim_value1{t},State_vars); 
                    
        %Probability of authoritarian parenting;    
        PS_prob = (1./(1+1./exp(value1-value0)));            
        PS = NaN(size(value1));
        %Simulate the behavior given random numbers ( fp_parfor.PS_forward );
        PS(isnan(PS_prob)==0) = ( PS_prob(isnan(PS_prob)==0)>= fp_parfor.PS_forward(isnan(PS_prob)==0,1,t) ) ;         
        
        %Optimal parental investments by parenting style (P=0);       
        inv0  = Interp_function_simulation(fp_parfor.estim_policy0{t},State_vars);
        inv0(inv0<fp.Min_Time)=fp.Min_Time;
        inv0(inv0>fp.Max_Time)=fp.Max_Time;

        %Optimal parental investments by parenting style (P=1);               
        inv1  = Interp_function_simulation(fp_parfor.estim_policy1{t},State_vars);
        inv1(inv1<fp.Min_Time)=fp.Min_Time;
        inv1(inv1>fp.Max_Time)=fp.Max_Time;

        %Optimal parental investments;
        Inv = inv0.*(PS==0) + inv1.*(PS==1);

        %Evolution of a child's skills given optimal parental behavior;
        Skills(:,t+1) = technology(t,alpha, hc{t}, Inv, PS, H_peers,  TFP );    
        %Store distribution of skills;
        hc{t+1} =  Skills(:,t+1);
            
        %Updating the grid points of the model when solving backward;
        grid{t} = prctile( hc{t} , fp.percentiles_kid ) ;
        grid_peers{t} = prctile( Peers(:,t) , fp.percentiles_kid ) ;
        grid_inv0{t}   = prctile( inv0 , fp.percentiles_inv ) ;
        grid_inv1{t}  = prctile( inv1 , fp.percentiles_inv ) ;

        %Store the aggregate state variables (distribution of skills and
        %parenting style in the economy);
        distrib_skills(:,t+1) =  Skills(:,t+1) ;
        distrib_PS(:,t) = PS;                               
        

        %%%%%%%%%%%%%%%%%%%%%%%%%;
        %Store Simulated Data;
        %%%%%%%%%%%%%%%%%%%%%%%%%;

        ind_obs = 1+(t-1)*fp_parfor.n_simulation:t*fp_parfor.n_simulation;
    
        sim_data( ind_obs , 1 ) = 1:fp_parfor.n_simulation;
        sim_data( ind_obs , 2 ) = repmat(fp.ages(t),length(ind_obs),1);
        sim_data( ind_obs , 3 ) = h;
        sim_data( ind_obs , 4 ) = hc{t};    
        sim_data( ind_obs , 5 ) = hc{t+1};    
        sim_data( ind_obs , 6 ) = Inv;    
        sim_data( ind_obs , 7)  = Peers(:,t) ; 
        sim_data( ind_obs , 8 ) = Peers_log(:,t);
        sim_data( ind_obs , 9 ) = PS; 
         
        sim_data_network( ind_obs , 1 ) = 1:fp_parfor.n_simulation;
        sim_data_network( ind_obs , 2 ) = repmat(fp.ages(t),length(ind_obs),1);
        sim_data_network( ind_obs , 3 ) = h;            
        sim_data_network( ind_obs , 4 ) = hc{t} ;  
        sim_data_network( ind_obs , 5 ) = n_friends;

    
end %t

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Add next-period peers into the simulated dataset;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
for t=1:(fp.periods-2)

        ind_obs = 1+(t-1)*fp_parfor.n_simulation:t*fp_parfor.n_simulation;
        sim_data( ind_obs , 10 ) = Peers(:,t+1); 
        sim_data( ind_obs , 11 ) = Peers_log(:,t+1); 

end

        sim_data( : , 12 ) =  fp_parfor.j ;


sf.distrib_skills = distrib_skills;
sf.distrib_PS = distrib_PS;
sf.grid = grid; 
sf.grid_peers = grid_peers;
sf.grid_inv0 = grid_inv0;
sf.grid_inv1 = grid_inv1;
sf.sim_data = sim_data;
sf.sim_data_network = sim_data_network;

